SEPARATION OF VARIABLES IN PERTURBED 

CYLINDERS 

A. ASLANYAN AND E.B. DAVIES 

Abstract. We study the Laplace operator subject to Dirichlet 
boundary conditions in a two-dimensional domain that is one-to- 
one mapped onto a cylinder (rectangle or infinite strip). As a result 
of this transformation the original eigenvalue problem is reduced 
to an equivalent problem for an operator with variable coefficients. 
Taking advantage of the simple geometry we separate variables by 
means of the Fourier decomposition method. The ODE system 
obtained in this way is then solved numerically yielding the eigen- 
values of the operator. The same approach allows us to find com- 
plex resonances arising in some non-compact domains. We discuss 
numerical examples related to quantum waveguide problems. 



1. Introduction 

The object of this study is the Dirichlet Laplacian in a perturbed 
cylinder, i. e. a domain that is mapped onto a rectangle or an infinitely 
long strip depending on the domain being compact or non-compact. A 
typical example of a perturbed cylinder is a waveguide where the propa- 
gation of waves is governed by the Helmholtz equation. The two major 
types of waves observed in waveguides are referred to as trapped modes 
and resonance solutions. Mathematically trapped modes are described 
as the eigenfunctions of self-adjoint operators. These are C? functions 
associated with bound states, as opposed to resonance solutions that 
correspond to the so-called scattering poles, or complex resonances. 
Despite their different nature, eigenvalues and resonances are some- 
times closely connected with each other. For example, eigenvalues may 
generate resonances under small perturbations of the domain. It is this 
situation that interests us and motivates our study of resonances. 
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Waveguide phenomena are usually associated with either Dirichlet or 
Neumann boundary conditions. The former correspond to scattering 
problems in quantum theory, the latter appear in acoustics. We refer 



to the series of papers ]T3|, [11], [L2| where bound states and scattering 
in quantum waveguides are studied. In |T7) the authors address the 
problem of finding quantum resonances numerically for a particular 
waveguide. For results on acoustic waveguides, both theoretical and 
numerical see, for instance, [JRI [TJL 0, 0, in KB. The papers cited 



here are concerned with either eigenvalues or resonances occurring in 
waveguides under specific conditions, or both of these. Our paper is 
very close in spirit to || and In the former the main issue is the 
resonance-eigenvalue connection, and the technique of the latter is also 
based on the separation of variables. 

Our intention is to study the above mentioned problems numerically. 
Dealing with both of them involves solving a boundary value problem 
for the Dirichlet Laplacian in two dimensions, that is either a self- 
adjoint eigenvalue problem or a non-self-adjoint resonance problem. 
Along with general numerical methods applicable in two dimensions, 
there exist techniques especially designed for cylinder-like domains also 
called ducts in acoustics. We have already mentioned |], [L5| where such 
methods are developed. Both of these papers stress the importance 
of advanced methods specially designed for acoustic waveguides. It is 
hardly surpsising that carefully performed numerical analysis is equally 
important for quantum problems. In [15] the authors apply a second- 
order finite difference method and implement an iterative procedure 
for the resulting algebraic system. It is mentioned there that standard 
methods not using any preconditioning are likely to fail especially when 
a large wave number is involved. 

The numerical approach proposed in is similar in spirit to that 
of our paper. In both cases the Helmholtz equation is reduced to the 
so-called coupled mode system of equations via the separation of vari- 
ables. The main difference is that in [|IJ the coefficients of the ODE 
system have to be computed numerically whereas our choice of the 
Fourier expansion functions allows us to find them in closed form. The 
Dirichlet problem studied here is separable as opposed to the much 
less straightforward Robin case. The examples in the cited paper are 
related to higher frequencies while we concentrate on the lowest oscil- 
lation mode only. On the other hand, the transfer method we use for 
the final ODE problem is able to handle a waveguide with a narrow 
throat — a situation not covered in [I]. 

The aim of this paper is to elaborate a method suitable for per- 
turbed cylinders that takes account of their geometry. The method 
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based on the Fourier decomposition in one direction allows us to sepa- 
rate variables in the Helmholtz equation explicitly leading to a system 
of ODEs. This is done for a fairly general geometry in the next sec- 
tion. In section 3 we discuss different boundary conditions involved. 
First we deal with standard self-adjoint conditions, then concentrate 
on non-compact domains and define resonances by a specific boundary 
condition at infinity. What is often called the radiation condition in the 
literature is rewritten in terms of the resulting ODE system. We end 
up with a non-self-adjoint eigenvalue problem on a finite interval whose 
solution approximates that of the original resonance problem. Finally, 
we use the transfer method of to find the eigenvalues of the two 
problems. Our numerical results illustrate the closeness of eigenvalues 
and resonances and are presented in section 4 where we conclude by 
discussing the rate of convergence. 

To be able to compare the method of the paper with others we look at 



the finite volume method described, e.g. in fig . As discussed in section 
H, the proposed approach tested on our eigenvalue examples proves 
to be significantly more efficient than the standard two-dimensional 
procedure. 

2. Separation of variables for the Laplacian 

2.1. Change of variables. Consider the operator H := —A acting 
on jC 2 (Q) subject to Dirichlet boundary conditions. The domain Q is 
defined as 

n = {{£,ri): a<£<b, < V < (2.1) 

in the Cartesian coordinates (£,77). The possibility of a and b being 
infinite is not excluded here, so that Q is not necessarily compact. The 
function </?(£) is assumed to be smooth and satisfy </?(£) > 0, £ G [a,b]. 
To find the spectrum of H we solve the Helmholtz equation 

-Af(t,r)) = \f(Z, V ), fori) en (2.2) 

with the boundary conditions 

= 0, (£,??) €<9fi (2.3) 

The change of variables 

x = £, y = ryMO (2.4) 

maps the perturbed cylinder Q onto = {(x,y) : a < x < b, < 
y < 1} which is either a rectangle or a strip (infinite or semi-infinite). 
We mention that a similar method has been used by Borisov et al. 
to study bound states associated with a local perturbation of a strip 
or layer. In our case the deformation reduces the width of the strip 
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locally, and there are no bound states. The transformation (|2.4J) can 
be expressed in the differential form as 



wv xy , v xy = ( f , w = ( n ^ 

0y / 







The quadratic form corresponding to H is given by 

j(f) = [ [(Ve,/, Ve,/) - A(/, /)] dedTy 

or, equivalent ly, by 

./(/) = / [(^V xy /, WV^/) - A(/, /)] <p(x)dxdy. 



'Ho 

This can be rewritten as 



J(/)= / [(^V^/, V,,/) - A(/, /)] ^(x)dxdy, (2.5) 



where 



A = W*W= [ _^ (1+ J{ )2) | • (2.6) 

Hence the Helmholtz equation in the new variables takes the form 

{iff x ) x - (ip'yf x ) y - Wyf v ) a + ((l + (<p'y) 2 )f y /ip) y + Xipf = 0. 

(2.7) 

Since we restrict ourselves to the Dirichlet case no change in the 
boundary conditions is required here, the condition (|2.3| ) is retained 
on cftV However, in a generic situation one can still use ( |2.4|) , ( |2.7| ) 
provided obvious changes are made to the original boundary condi- 
tions where necessary. For instance, instead of Neumann boundary 
conditions at rj = ip(£) one would have 

while Neumann boundary conditions at x = a, b would become 

(<pf x -<p'yf y )\ x=a)b 

Note that a similar transformation can be also done for a more gen- 
eral domain 

ft = {(£,7?) : a < £ < b, <^(0 < 77 < <^ 2 (0}- 
The change of variables 

c V-<Pi 
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leads to a quadratic form of type ( |2.5| ) whose coefficients are not quoted 
here for brevity. 



2.2. Discretisation in the y-direction. The quadratic form ( |2.5| ) is 
related to the transformed operator on the weighted space C 2 (Q , <pdx dy). 
Here and below we use the notation V = V ' xy . 

To discretise the form J(f) in the ^/-direction let us separate the 
variables expanding / as 



oo 

£ 

k=l 



gk{y)h k {x) 



(2.8) 



Recall that we have Dirichlet boundary conditions everywhere so that 
our natural choice is to work with an orthonormal system of functions 
vanishing at the horizontal parts of the boundary. We therefore opt for 

g k = \/2sin(7r%) . 

Denote f = f, fi = f x , f 2 = f y ; 



k=l 



where 

9° k = 9k = 9k, g\ = v / 2cos(vr%). 
In this notation the x-dependence is determined by the functions 



h° k 



hi 



hi = hi 



hi 



Trkhi 



L k "'ki 

(throughout the paper ' denotes differentiation with respect to x). 
We notice that the variables are separated in the coefficients of </(/): 

vWAifay) = B i3 (x) + C l] (y)D ij (x), i,j = 1,2. (2.9) 

The entries Ay of the matrix A are defined by ( |2.6| ); the matrices B, C 
and D satisfying the above decomposition are given below: 





(f 
i 



B = 
The formula 



i 



C 



y 



D 



J(f) 



y y j v -v 5 

allows us to rewrite (|2.5|) in the form 

' 2 



£_ 

f 



2 



dxdy 



fifjdy + Dij / Cijfifjdy-Xip \f\ 2 dy 



dx 
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BijEij + DijFij — \{pG dx. 

As we substitute the expansion (|2.8|) into the above integral, the coef- 
ficients Eij, Fij, G are readily computed below. 

Jo Jo 



k,r ^ k,r 



i.l 



Cijfifjdy = [ dj V hig* k V Kg^dy 
Jo k 

= f CMdy = 52$MK'> 

k,r ^° k,r 

JO t, „ i. 



k r k 

To find a^ r we use the orthogonality relations for g % k \ in fact, we only 

33 
X kr 



need the diagonal elements a kr = 6 kr - The coefficients /3 kr are also 



calculated in the closed form: 



Pkr ~ u ' ykr 



I j 1 

3 ~ 2vr 2 fc 2 ' 

k+r 4(fc 2 +r 2 ) 



-1 



= r 

7r 2 (fc2 _ f ,2)2, 7^ r 



nl2 _ o^l 
r^fcr rVfc 



21 



k = r 



l 

'2?rfc' 

\fc+r 2A; r, / 

r r 2_ fc 2p «- 7^ ' 



7r(r 2 — fc 2 ) : 

The quadratic form is now reduced to that of a one-dimensional differ- 
ential problem. 

2.3. Canonical ODE system. Having done the above calculations 
we finally arrive at 



J(f) 



+ Y ^ 2 kr^h k h r - ^vr (rpl 2 r h' k h r + kfi? k h k K r ) 



E i^i 2 + 



Tik 



99 dx. 
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The Euler equations are easily derived in the standard way. A simple 
calculation shows that the discretised form ( |2.3|) is equivalent to the 
ODE system written in its canonical self-adjoint form as 



-{Ph')' + Qti - {Q*h)' + Rh = 0. 
Here the vector of unknowns 

h = (h 1 ,h 2 ,...) T ; 
the matrix coefficients are given by 

Pkr = <fhr, Qkr = -irkPllip' 



(2.10) 



R, 



kr 



n 2 kr/3 f 



22 f 
kr 



12 



Jkri 



k,r= 1,2, 



For practical purposes we truncate the system to a finite number of 
equations taking a sufficiently large N and keeping the same notation 
h, P, Q, R for the truncated matrices where k,r = 1, . . . , N. This is 
justified by the fact that the Fourier coefficients involved in ( |2.8|) are 
rapidly decaying in k and therefore higher order terms can be neglected. 
In |] it has been suggested that N should be of order h\/X where h 
denotes the mean width of the duct if the curvature of its boundary is 
not too large. In the examples of section f| the width of the waveguide 
varies greatly from point to point. The size of iV is determined exper- 
imentally and found to depend mainly on the width of the narrowest 
portion of the waveguide. 

Equivalently, we reduce ( |2.10| ) to the Hamiltonian system of 2N 
equations: 



JH' = K(x, X)H, x G 



(2.11) 



where 



H 



h 

Ph' + Q*h 



e C 1 



J 



o -/ 

/ 



K 



-R + QP- l Q* 
-P~ l Q* 



QP- 
P- 1 



The system (pip) is self-adjoint with P = P* > 0, R = R* for A G R. 
We can therefore apply advanced numerical methods (see, for example, 
0, 0|) to find the eigenvalues of the problem and the relevant solutions. 
Before proceeding to this task let us discuss the issue of boundary 
conditions. 



s 
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3. Boundary conditions 

3.1. Self-adjoint problem. To make sure the original boundary con- 
ditions are involved in the ODE problem consider a generic situation 
when we have a functional 

J{f)= F(x,y,f,f x ,f y )dxdy. 

To derive the corresponding Euler equation we replace / by / + £7 and 
compute 

+ e I 7 {F u dy - F fy dx) . 

JdQo 

Putting the first integral equal to zero we obtain the differential equa- 
tion ( |2.7|) ; the second is responsible for boundary conditions. For our 
class of problems dQo = {y = 0} U {y = 1} U {x = a} U {x = b}. The 
conditions at different parts of the boundary are defined by 



J a 

[\F f J dy = 0. 
Jo 



Taking account of the obtained quadratic form, we get 

lfy\ y =o dx = °> (3- 1 ) 



a 



7 (tf(x)fx ~ 1 + <f ft 
V <P 



dx = (3.2) 

y=l 



on the horizontal lines. As pointed out in subsection 2T, Dirichlet 
boundary conditions remain unchanged in the new variables and are 
automatically taken into account by virtue of our choice of the functions 
gk(y) in (|2.8| ). The above integrals (|3.1|) , (|3.2|) vanish because of the 
implied condition 7 = 0. A difficulty would only occur if we had 
more complicated conditions at the curvilinear part of the boundary 
of Q. Dirichlet boundary conditions are the ones relevant for quantum 
mechanical problems and they enable us to separate the variables in the 
quadratic form explicitly. We refer to [1]] where the authors consider 
arbitrary boundary conditions of the form (af + b^j \ QQ = by using 
appropriate orthogonal curvilinear coordinates. The problem of this 
kind requires a more complicated expansion to be used instead of 
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In that case Fourier coefficients are not obtained in closed form but 
should be calculated numerically. 

On the vertical parts of the boundary we have 



/ 9 (¥>fx-v'yf y )\ x=aib dy = o. 

JO 



(3.3) 



The Dirichlet case is as easy to treat as before: the conditions hk{a) 
hk{b) = 0, k = 1, . . . , N are imposed on the solutions of ( 2.11j) . Con- 



sider also a domain where y?'( a ) — v'ify = — the situation typical 
for compactly perturbed strips and, in particular, for some waveguides. 
Here we are able to handle a more general case. For instance, Neu- 
mann boundary conditions at x — a, b do not change and become 
h' k (a) = h' k (b), k = 1, . . . , N in terms of the system (|2.11|). However, 



one cannot fully separate variables in generic Robin conditions of form 



3.2. Radiation condition. Bearing in mind the resonance problem 
that is of main interest to us, let us consider a domain Q such that, in 
notation of subsection 12711 x G (a, oo) and 



<p(x) ~ 1, x > X (3.4) 

for some X > 0. Similar assumptions are often made in papers dealing 
with scattering problems, for instance in [l], |15| ||. We take Dirichlet 
boundary conditions on dQ and require a different type of condition to 
be satisfied as \x\ — > oo. Namely, for a given A G C there always exists 
a unique solution of ( |2.2|) that has the form 



k=2 



f(x, y) = (exp(-t!x) + st exp(tirr)) g x {y) + ^ s k9k{y) exp(t k x), x > X. 

(3.5) 



t k = -^(nk) 2 - A, Ret fc < 0; 



Here we denote 



g k are the same as in subsection |2.2j . The coefficients s k , k = 1,2, 
are defined by the formula (|3.5| ) uniquely for each value of A. We put 
uj = a/A and consider s k as functions of u. The function si called 
the scattering coefficient of the problem is involved in the definition of 
resonances. The reader will find their general definition in ||. Note 
that when variables are separated the following construction proves to 
be more handy. If the scattering coefficient s\{uj) has a pole at u = u>o 
we say that A = Uq is a resonance. This does not include all the 
resonances but only those lying on the first non-physical sheet (see 
for detailed explanation). We refer to |T7|, Q for the equivalence of the 
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two definitions. Apart from its simplicity, the approach based on ( |3.5| ) 
has another distinctive advantage. It is known from scattering theory 
that Si(uj)si(uj) = 1 and that Si is analytic in the half-plane Imu < 0. 
Therefore instead of seeking the poles of si(u>) one can look for its 
zeros located in the lower tu-half-plane. This is the approach we use 
here along with the separation of variables in the perturbed cylinder 
tt. 

Given the above definition, there is an obvious difference between the 
resonance problem and a classical spectral problem. Indeed, according 
to ( |3.5p here we are looking for a solution exponentially growing at 
infinity. Note, however, that the resonances we are interested in occur 
as perturbations of eigenvalues and are typically situated near the real 
axis. This means that the values of |Imifc|, k = 1,2,..., are rather 
small and therefore the corresponding solution grows slowly. 

Combining ( |3.5|) with (|2.8| ) for a sufficiently large x we get 

h\(x) = exp(-tx) + siexp(tx), h' k (x) = 4/^(2), k — 2, 3, ... . 

There are two ways to handle these conditions. One can solve the 
inhomogeneous problem (as has been done in |5| in two dimensions), 
then find the zeros of S\. Alternatively one can put S\ = straight 
away, then solve the resulting eigenvalue problem with A-dependent 
boundary conditions. The latter approach leads to the set of boundary 
conditions at X 

i>xH(X) = 0, i) X = (T-P- l Q\ P- 1 ), T = diag(t!,-t 2 ... ,-t N ). 

(3.6) 

This formula known as the radiation, or outgoing wave condition singles 
out the solution whose first component grows and the others decay 
exponentially at infinity. It is this solution that is sometimes called the 
resonance eigenfunction. 

It should be mentioned here that our approach agrees with that using 
exterior complex scaling (see, for example, ||). In this technique one 
replaces the operator by a family of operators on the same domain 
which depend analytically on a complex parameter. The operators are 
independent of the parameter for x < X and are associated with a 
space scaling for x > X. One computes the complex eigenvalues of 
this family acting in £ 2 (Q) and proves that they do not depend on the 
parameter, subject to certain conditions. It is known that the complex 
eigenvalues coincide with resonances defined via either the scattering 
coefficients or analytic continuation of the resolvent kernels. One may 
verify that exterior complex scaling yields the same boundary condition 
at x = X as fl3lf). 
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4. Numerical examples 

4.1. Transfer method. Summarising the results of the first three sec- 
tions let us formulate the problems to be solved numerically. We are 
looking for such values of A that the system fl2.11|) has a non-trivial 
solution satisfying 

ij a H(a) = 0, i> b H{b) = 

where 

1 ^ a = A = (/, 0); 

2 ip a = (J, 0), ip b = ip x as defined by (U). 

Problem 1 provides approximations to the Dirichlet eigenvalues of a 
compact domain of kind ( |2.1| ); problem 2 enables us to calculate com- 



plex resonances that may occur in an unbounded domain of the same 
type satisfying fl3.4|). 

The method we apply to both problems is based on the orthogonal 
transfer of || that we shall briefly outline below. The manifold of the 
solutions of the system (|2.11|) satisfying the left boundary condition is 



determined by 

ip(x)H(x) = 0, a < x < b, 
where ip £ <C Nx2N solves the Cauchy problem 

■0' = ipJK, ip(a) = ip a . 
Theoretically one can integrate the above equation for a fixed A, define 



/(A) = det 



f^(b;X) 

V <MA) 



and solve /(A) = to find the eigenvalues of the problem. This method 
is known to be hopelessly inefficient because ip(x), although formally 
of rank N, can have almost linearly dependent rows. Abramov 
proposed replacing ip by ip(x) = v(x)ip(x), where v e C NxN , det v ^ 0. 
The function v is chosen so that to ensure ip{x)ij)*ix) = const. The 
transfer equation now takes the form 

$ + 4>JA (i - vi*(00TV) = °> ^(«) = Vv (4.1) 

The RHS of ( |4.1| ) is bounded and the solution tp exists on the whole 
of [a, b}. By comparison with ip, the matrix if; has the key advantage 
of being easily computed without loss of rank. The use of this idea 
proved essential to obtain stable results for this problem. 

Having calculated the smooth function ip we proceed to find the 
eigenvalues. For the resonance problem this is done along the lines of 
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H where the idea of || has been applied to non-self-adjoint eigenvalue 
problems. As observed there, 

/(A) = detf ^ ^ /(A )det^ 

so that the zeros of / and / coincide. Moreover, the zeros of / can be 
found by the method based on the argument principle although / does 
not have to be analytic in A as opposed to /. Still the number of zeros 
of / inside a contour T is 

N = i- j> dArgf(X) 

as shown in @]. It is this computational formula that we use to locate 
the complex eigenvalues of problem 2. Taking a shrinking sequence of 
contours T we find the zeros up to a chosen accuracy. We computed 
the contour integrals reliably for circles of radii down to 10 -4 and made 
sure that if the centre was shifted by a similar order of magnitude 
the integrals vanished. In the better conditioned case 1 we applied 
Newton's method allowing us to calculate the eigenvalues of the self- 
adjoint problem. 

Note that typically problem 2 is much harder to solve than 1, and 
our examples are no exception. When resonances are situated near the 
real axis the scattering coefficient s\ has a pole and a zero close to one 
another. Naturally, the closer they are the less stable the problem is. 

4.2. Results of computations. As an example we consider a quan- 
tum waveguide with indentations defined in the Cartesian coordinates 
(£,ti) as 

W = {-oo < £ < oo, 0<7]< <p(g) = 1 - a (V (? ~ 7)2 + e~ (?+7)2 ) } 

where a and 7 are real positive constants (see fig. 1.1). 
We shall be working in two different domains: 

n x = w n {0 < £ < 7}, o 2 = w n {o < £}. 

The domains fli and relate to problems 1 and 2 of the previous 
subsection, respectively. In the latter the resonance boundary condition 



is imposed at a sufficiently far point X as suggested in subsection p?2 . 
In our experiments we put 7 = 2, so that it suffices to take X > 5.4 to 
ensure <p'(X) < 1(T 4 . 

In this example there exists a = a* ~ 1 such that the two parts 
of dQ touch one another near £ = 7. For this value of a the do- 
main Q consists of three disjoint parts as shown in fig. 1.2, so that the 
eigenvalue problem is decomposed into three separate problems. The 
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Laplacian considered in the compact domain has infinitely many 
real eigenvalues accumulating at infinity As we decrease a joining the 
three subdomains, one expects the eigenvalues to disappear generat- 
ing resonances in their neighbourhood. A similar phenomenon where 
resonances originate from eigenvalues as the domain is perturbed has 
been observed in B] although the mechanism by which they emerge is 
different here. 

Clearly, both domains Q\ and Q2 satisfy the conditions of subsection 
PTT| for a < a*. We compute resonances, i.e. eigenvalues of problem 
2 for a range of a using the perturbed cylinder approach. In the self- 
adjoint example the eigenvalues of problem 1 are found for the same 
values of a. The ODE problem is solved by the transfer method de- 
scribed in the previous subsection. An important question is how to 
choose the number of terms to be retained in (2.8) or, in other words, 
the dimension of the system ( 2.1 1|) . This number should depend on a: 
indeed, for a close to a* the width of the domain is small near £ = 7 so 
that one needs to keep a larger number of terms Qk{y) in (|2.8|) . There 
are two possibilities here: first, to increase repeatedly the number of 
terms by one and solve the system with the corresponding constant 
number of unknowns until the answers converge. Secondly, instead of 
keeping a large number of terms throughout the interval one can start 



off with a smaller N in ( |2.11| ) . Moving along the x-interval one changes 
N gradually adding or removing variables hk(x) depending on the size 
of <p(x). We have been using both techniques in different situations 
ensuring that the results coincide within the chosen accuracy for two 
subsequent values of N. For reasonably small values of a it suffices to 
take smaller N: for instance, when a = 0.8 the results obtained for 
N = A and greater coincide up to the tolerance of 10~ 4 . The maximal 
number of terms taken in our computations is N = 30 for a = 0.97 
(the largest value considered). 

There is an important connection between the two problems 1 and 
2 which makes us study them within the same framework. Namely, as 
a — > a* both the eigenvalues of Q\ and the resonances of Q2 converge 
to the Dirichlet eigenvalues of the domain Q*. These cannot be found 
by the same method because fl* has a cusp at £ = 7 and is not a 
perturbed cylinder in our terminology. In this case one could still sep- 
arate variables in the Helmholtz equation in a similar way arriving at 
a singular ODE problem. This question requires a separate considera- 
tion which is beyond the subject of our paper. Recall that the aim of 
our numerical experiments is to calculate resonances occurring in this 
example and find out how they are related to eigenvalues. To be able 
to make proper comparisons we have used the finite volume method 
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to discretise the operator on and find its spectrum. To know the 
limit eigenvalue is also helpful as it serves as an initial guess for the 
eigenvalues of problems 1 and 2. 

We applied the finite volume method to a series of problems of type 
1 to compare the effectiveness of the two approaches. As expected, the 
comparisons are in favour of the discussed method, which appears to 
be several times faster than the conventional one. The benefits of our 
technique are more spectacular for larger values of a. For example, 
when a = 0.8 it takes twice as long to get accurate results by the finite 
volume method, while for a = 0.9 the new method is almost four times 
faster. Everything else being equal, the closer a gets to its critical value 
the more advantageous the perturbed cylinder approach is. 

The numerical results presented below are related to the lowest 
Dirichlet eigenvalue of f2* and quoted in terms of the wave number 
u = y/\ (we shall retain the term eigenvalue for the wave numbers). 
The smallest eigenvalue corresponding to the domain with the cusp is 
co>* = 4.6252; it is included in the diagrams to illustrate the convergence 
of our results. 

In fig. 2 a series of the eigenvalues of problem 1 for a € [0.7, 0.97] 
is shown. They converge to a;* as a — > a*. The linear rate of conver- 
gence is in agreement with standard perturbation theory: the principal 
correction term is of order e = a* — a since the coefficients of ( [2.1 1|) 
depend on a linearly. 

The resonances u; 2 emerging from the lowest eigenvalue as a de- 
creases can be found in table 1. They are also shown in fig. 3 where 
their real parts are plotted against their imaginary parts. They def- 
initely converge to cu*. The problem is, of course, very sensitive to 
perturbations near the critical value of a = a*, or e = 0. However, 
we believe that for a < 0.97 our computations are relatively stable 
and provide reliable results. For larger values of a we have observed 
various instability effects preventing us from computing resonances ac- 
curately. The system ( |2.11| ) is difficult to solve for values of uj close 
to oj* (and, consequently, the corresponding resonance) because its co- 
efficients change very rapidly for a ~ a*. Even methods suitable for 
stiff systems fail to produce satisfactory results when e is too small. As 
discussed in the previous subsection, another reason why the problem 
is likely to be unstable is the closeness of resonances to the real axis. 
These two factors make calculations slow and inefficient for a close to 
a*. 

The a-dependence of the real and imaginary parts of the resonance 
originating from u;* is shown in fig. 4. The real part is found to depend 
on the perturbation parameter linearly, whereas the imaginary part 
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Table 1 . Resonances L02 



15 



a 


UJ2 


a 


Ll>2 


0.7 


4.2988 + 0.0545i 


0.9 


4.5250 + 0.0050i 


0.75 


4.3715 + 0.0348i 


0.925 


4.5492 + 0.0032i 


0.8 


4.4223 + 0.0212i 


0.95 


4.5755 + 0.0017i 


0.825 


4.4498 + 0.0154i 


0.96 


4.5864 + 0.0012i 


0.85 


4.4741 +0.0113* 


0.97 


4.5967 + 0.0008i 



Table 2 . The values of 5 = u\ — Re 002 



e 


0.3 


0.25 


0.2 


0.175 


0.15 


5 


0.0436 


0.0256 


0.0182 


0.0144 


0.0123 


e 


0.1 


0.075 


0.05 


0.04 


0.03 


5 


0.0081 


0.0044 


0.0015 


0.0008 


0.0004 



seems to be of order e p with p « 3/2 for e < 0.2. A different behaviour 
has been observed in ||, where for a different geometry the authors 
derived an asymptotic formula for the imaginary part of the resonance. 
The resonances in the example of that paper are shown to be analytic in 
the domain perturbation parameter, their imaginary parts depending 
quadratically on the parameter. 

One generally expects that if one perturbs an eigenvalue which is 
embedded in the continuous spectrum, then it is transformed into a 
resonance near the real axis. In many cases one can even write out a 
perturbation expansion, and a general theory for some such cases was 
described by Agmon . However in our situation the natural parame- 
ter e can only take positive values for obvious reasons. Therefore even 
if the resonance depends analytically on e for e > 0, there is no proof 
that it has an expansion with finite coefficients around e = 0, nor even 
that the resonance converges to the eigenvalue as e j 0. The numerical 
experiments do, however, suggest that not only does it converge, but 
also that its real part and imaginary parts are expanded in powers of 

In the table below the eigenvalues of problem 1 are compared with 
the real parts of the resonances related to the same values of a. We 
tabulate the difference 6 = u)\ — Re 102 between the Dirichlet eigenvalue 
and the real part of the associated resonance and observe that this 
quantity decreases as e | 0. One might hope to deduce the rate of con- 
vergence of 5(e) from asymptotic perturbation formulae. The question 
how to obtain such formulae for the resonance still remains open. 
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5. Conclusions 

The main idea of this paper was to reduce an eigenvalue problem in 
a two-dimensional perturbed cylinder to an ODE problem. Complex 
resonances occurring in perturbed strips were also dealt with in the 
same way. 

The method of this paper allowed us to discretise the problem in 
one direction taking into account the geometry of the domain. This 
was especially important for the irregularly shaped domains of section 
U Standard finite difference methods would require a significant mesh 
refinement in the narrow part of the considered waveguide. For com- 
parison purposes we also computed the eigenvalues of problem 1 by the 
finite volume method. As discussed in the previous section, it proved 
to be substantially more time-consuming than the method based on the 
separation of variables for the class of problems studied here. This is 
in agreement with the already mentioned results of [|], [15| where other 
advanced methods were designed to suit similar problems. We believe 
that our approach is competitive and recommend it for ill-conditioned 
eigenvalue and resonance problems. Our confidence is supported by 
the strong agreement between numerical results and analytic expecta- 
tions. Indeed, the convergence of both eigenvalues and resonances to 
their limit value (computed by an independent method) confirms the 
reliability of the proposed technique. 
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1) a = 0.9 




2) a = 0.9999 
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Figure 1. Waveguide with indentations 




Figure 2. Eigenvalues of problem 1 for a range of a 
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Figure 3. Resonances generated by u;* 
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Figure 4. Real and imaginary parts of resonances: a-dependence 
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